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f*^ 1 Abstract 

In this note we discuss the construction of high order asymptotic preserving numerical schemes for the 
, Boltzmann equation. The methods are based on the use of Implicit-Explicit (IMEX) Runge-Kutta 

■ methods combined with a penalization technique recently introduced in [6j. 
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Resume 

Dans cette note nous discutons la construction de schemas d'ordre eleve pour l'equation de Boltzmann 
qui preservent la limite asymptotique. Les methodes sont basees sur l'utilisation de schemas de Runge- 
Kutta explicites-implicites combinees avec une technique de penalisation introduit recemment par [BJ. 

Mots-cles : Methodes Runge-Kutta Implicites-Explicites, equations raides, equation de Boltzmann, 
X/ , limite fluide, schemas preservant l'asymptotique. 
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Les equations cintiques, comme l'equation de Boltzmann sont utilise avec succes dans de nombreuses 
applications reelles. L'equation de Boltzmann decrit revolution temporelle de la fonction de distribu- 
tion d'un gaz avec des interactions binaires elastiques. II est importantc de mentionner que la solution 
numerique de l'operateur de collision represente un defi majeur pour les methodes numeriques tradition- 
nelles qui n'est pas encore resolu. Cela est particulierement vrai en proximite des regimes fluidcs. Dans 
ces regimes le taux des collisions intermoleculaires crot de faon exponentielle et done le temps entre deux 
collisions successives devient tres petit. D'autre part, l'echelle de temps reel pour revolution du gaz est 
l'echelle de temps de la dynamique des fluides, qui est normalement beaucoup plus grande que le temps 
entre deux collisions. Une mesure de l'importance des collisions est donnee par le nombre de Knudsen 
e, qui est grand dans la limite rarefiee et petit dans la limite fluide. Ainsi, les approches numeriques 
traditionnelles perdent leur emcacite en raison de la necessite d'utiliser de temps tres petits pour la dis- 
cretization temporelle. Nous rappelons que la discretisation directe en temps de l'equation de Boltzmann 
est un gros probleme dans les regimes raides en raison de la haute dimensionnalite et de la non-linearite 
de l'operateur de collision qui rend peu pratique l'utilisation de solveurs implicites. 

Plusieurs auteurs ont aborde le probleme dans le recent passe( voir [H [SJ |H1 HI M HI] ) et les references 
a l'interieur). Une strategie, parmi les plus puissantes, consiste en la construction des schemas dits 
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preservant l'asymptotique. Ces techniques permettent de resoudre le probleme dans tout le domaine pour 
tous les choix de pas de temps et de nombre de Knudsen. Dans cette note, nous proposons une nouvelle 
classe de schemas Runge-Kutta Implicites-Explicites pour l'equation de Boltzmann. Pour construire nos 
schemas, nous considerons une decomposition du terme de gain de l'operateur de collision en un partie 
en equilibre et en une partie en non equilibre. Cette decomposition de l'integrale de Boltzmann a ete 
egalement introduite par Jin et Filbet dans [B]. Les principaux avantages de l'approche proposee ici 
est que cela fonctionne de maniere uniforme pour une large gamme de nombres de Knudsen et evite la 
solution d'un systeme d'equations non lineaires, meme dans les regimes raides. De meme que pour [4], 
nous obtenons des conditions suffisantes pour la stabilite asymptotique et la preservation asymptotique 
de l'ordre temporel des schemas. En plus, nous construirons les schemas tels qu'ils preservent la positivite 
des solutions et les quantites physiques conservees. Pour plus de details nous renvoyons a [5]. 



1 Introduction 

The computation of fluid-kinetic interfaces and asymptotic behaviors involves multiple scales where most 
numerical methods lose their efficiency because they are forced to operate on a very short time scale 
(see [H El El El El E] and the references therein for a more complete bibliography). The Boltzmann 
equation close to fluid regimes represents the prototype example 

d t f + v-V x f=-Q(f,f). (1) 

£ 

Here f{x,v,t) is a non negative function describing the time evolution of the distribution of particles 
with velocity »el 3 and position x e Q C M. d * at time t > 0. 

The operator Q(f, /) describes the particles interactions. In the general case of the Boltzmann binary 
collision, it has the form 

Qb (/,/)= / B(\v-v m \,nMvV(vl)-f(v)f(v.)]dv m dn (2) 

where 

v' = v + -(v- «») + -\v - v*\n, v'„ = v + -{v - u») - - \v - u*|n, (3) 

and B(\v — v*|,n) is a nonnegative collision kernel characterizing the details of the collision. 

The Knudsen number e > is a non dimensional measure of the importance of collisions and is 
large in rarefied regions and small where the system is close to the fluid limit. In the latter regime, the 
intcrmolecular collision rate grows quickly and thus the collisional time scale becomes very small. On the 
other hand, the actual time scale of the evolution is the fluid dynamic scale, which can be much larger 
than the collisional time. 

In fact, for small values of e the distribution function is well approximated by a local Maxwellian 

M [f) = —P e ^(±^L], (4) 
U 1 (27tT) 3 / 2 v \ 2T J ' K ' 

where p, w, T are the density, mean velocity and temperature of the gas in the x-position and at time t 
defined as 

(p,pw,E) T = J ffl, v ,^\ dv, T = - p\w\ 2 ). (5) 

Now, passing to the limit for e — > and integrating ([T]) against 1, v and v 2 we recover the system of 
compressible Euler equations 

dtu + V x ■ F(u) = (6) 



2 



with 

u=(p,w,E) T , F(u) = (pw, gw <£> (w + pi), Ew + pw) T , p = pT, 

where I is the identity matrix. 

Implicit-Explicit (IMEX) Runge-Kutta schemes represent a powerful tool for the numerical treatment 
of stiff terms in PDEs |H |31 E]. When necessary they can be designed in order to achieve suitable 
asymptotic preserving {AP) properties. Their direct application to the Boltzmann equation however 
is not trivial since the complicated nonlinear structure of the collisional operator makes prohibitively 
expensive the use of implicit solvers for the stiff collision term. Additional difficulties are given by 
the need to preserve the most relevant physical properties of the solution, like conservation of mass, 
momentum and energy, nonnegativity and entropy inequality. In this short note we will illustrate how 
the introduction of a suitable penalization technique as in [5] permits to extend succesfully the IMEX 
formalism also to the challenging case of the Boltzmann equation. 



2 IMEX schemes for the Boltzmann equation 

In order to apply efficiently the IMEX Runge-Kutta approach to the Boltzmann equation we must avoid 
the prohibitive cost of the implicit evaluation of the stiff collision term. In order to achieve this we first 
reformulate the collision part using a suitable penalization term. 

2.1 Decomposition of the collision integral 

First, we observe that we can rewrite Qb(I, f) as 7 

Qs(/,/) = -(^(/,/)-M/), (7) 

e 

where P(f, /) = /) + pf and p > is a constant such that P(/, /) > 0. 

Observe that, by construction, the following property is verified by the operator P(f, f) 

Thus, P(f, f)/p is a density function and we can consider the following decomposition 

P(fJ)/H = M[f]+g, (9) 

where the function g represents the deviations from equilibrium of P( f, /). 
Thus the collision operator can be rewritten in the form 

QbU, f) = -9 + -{M[f\ -/) = - (W^- - M[f]) + »{M[f) - /). (10) 
e e e \ p, J e 

The above reformulation is equivalent to the penalization method for the collision operator recently 
introduced in [6j. Clearly, since the problem is stiff as a whole a fully implicit method should be used 
in the numerical integration to avoid stability constraints of the type At = 0{e). On the other hand, 
the linear part itself (M[f] — /) suffices to characterize the correct large time behavior of /. Therefore, 
instead of fully implicit methods, one may use methods which are implicit in the linear part and explicit 
in the non-linear part. This however, as we will see, introduces some additional stability requirements in 
order for the IMEX schemes to preserve the asymptotic behavior of the equation. 
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2.2 Application to the Boltzmann equation 

We can now introduce the general class of IMEX Runge-Kutta schemes for the Boltzmann equation in 
the form 

i-l 



F« = f n + At^H (£g(F&) - v ■ V*F«) + At^a^(M[F«] - F«) (11) 

3=1 3=1 

/ n+1 = /" + Ai^w, - « • V^FW) + AtJ2^{M[F^} - F«). (12) 



In the above scheme the explicit method is characterized by the i> x v matrix A — (Sy), a,j =0, J > i 
and the coefficient vectors are c = (ci, . . . ,c u ) T , Ci = ^ = (^li ■ • • j ^v) T > whereas the implicit 

method is a diagonally implicit Runge-Kutta (DIRK) defined by the v x v matrix A = {flij), aij = 0, 
j > i, and the coefficient vectors are c = (ci, . . . , c„) r , Cj = Sj=i a «i' w = (^l) • • - > w v) T ■ 
Let us first recall the definition of asymptotic preserving property [5] 

Definition 1 The IMEX scheme lllMty) for the Boltzmann equation is asymptotic preserving (AP) if in 
the limit e — > the scheme becomes a consistent discretization of the limit system of the Euler equations 

®- 

Note that if we multiply the IMEX scheme by the vector of collision invariants </>(v) = (1, v, v 2 /2) T 
and integrate in v we get a moment scheme characterized by the explicit method 

/ F®<f>(v) dv= f f n 4>{v) dv - At V aij [ v ■ V x F^<j>(v) dv (13) 
Jr 3 Jr 3 j=1 Jr 3 

f f n+1 (j)(v)dv = f f n <t>{v) dv - AiV^i / v ■ V x F (i U(v) dv. (14) 
Jr 3 Jr 3 i=1 Jr 3 

Thus a sufficient condition for a scheme to satisfy the ^4P property is that as e — >• we get F^ — > M[F^], 
i = 1, . . . ,v in (jlip . In addition we must require the kinetic numerical solution f n+1 to satisfy some 
additional numerical stability requirement. We illustrate this aspect in the sequel. 
First let us start with the following Lemma [5] 

Lemma 1 If all diagonal element of the triangular coefficient matrix A that characterize the DIRK 
scheme in equations are non zero, then 

limF (i) = M[F (i) l. (15) 

Formally Lemma 1 guarantees the AP property of the scheme. However, as opposite to the case of 
hyperbolic systems with relaxation, now because of the decomposition of the collision operator the last 
level (fi"2"l) still depends on e. After some manipulations it reads 



! n+l = r fl >_>,//,,, j - At^ffii • V B F« - \g{F^ 
+ AtJ2 Wihfijh (v ■ V X F^ - - g (F^)) + Y, WibijF®, 



(16) 



where bij are the elements of A^ 1 . The above expression turns out to be unbounded as e — > thus 
originating an unstable scheme. 

We introduce the following definition [2 [5] 
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Definition 2 An IMEX scheme in the form m\ )-kl 6 2 \) is globally stiffly accurate if the following condi- 
tions are satisfied 

Wi = a ui , Wi = a vi , Vi = l,...,v. (17) 

We can finally state the main result [5] 

Theorem 1 7/det A ^ and the IMEX scheme fllD - H^) is globally stiffly accurate, in the limit e — > 0, 
the IMEX scheme becomes the explicit RK scheme characterized by (A,w,c) applied to the limit Euler 
system 

In order to prove the Theorem it is enough to observe that the stiffly accurate property implies immedi- 
ately that f n+1 = fM. 

Remark 1 

• Theorem above guarantees not only asymptotic preservation but also asymptotic accuracy, namely 
the order of the scheme is preserved in the e — > limit. 

• The previous results can be extended to the case of CK-type schemes J3jj with an — 0, like the ones 
considered in Jffjj. However in this case asymptotic accuracy holds true only if the initial data are 
an 0(e) perturbation of the local Maxwellian equilibrium. 



2.3 Convexity of the schemes 

The determination of general conditions for positivity of the numerical solution in the space non homo- 
geneous case is quite difficult. Here we focus on the space homogeneous situation. Not that even in this 
case due to the reformulation of the collision term the analysis involve the whole IMEX scheme. Moreover 
the analysis here depends on the particular operator used as a penalization. 

Using the fact that in the space homogenous situation M[f] does not depend on time the IMEX 
scheme ([TT1) - (jl2D can be rewritten as 

A/" + 2 1 " + M \n fa ~ c h ) \ (18) 

j=i M J 

/( .«> . r + !^±^^^} + !^± Wl(Mlr] . F( ,y ) (19) 

£ * — ' a e * — ' 

i— 1 i— 1 

where A = ej (/iAt) and bij are the elements of (XI + A)^ 1 . 

The following theorem gives sufficient conditions for the above expression to represent a convex com- 
bination of probability densities. 

Theorem 2 A sufficient condition to guarantee that f n+1 > when f n > in Iil8]) -il9 \) is that the 
scheme is globally stiffly accurate and the following conditions holds true 

i i 

0<Y,*bihC h <l, Q<Y,bih{c h -c h )<l, Vi = l,...,!A (20) 

h=l h=l 

i 

0< ^2 b ih a hj <1, Vi = l,. j = l,...,i-l. (21) 
h=j+i 

Since the result is based on a convexity argument we also have an entropic result for the schemes. 
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Figure 1: L\ error for the distribution function / for the second (left) and the third (right) IMEX-BE 
method. 



Corollary 1 Under the assumptions of Theorem 1, if in addition the operator P(f, f) satisfies 

fp(fjy 



H 



<H(f), H(f)= f f log fdv, 
Jm 3 



(22) 



then H{f n+1 ) < H{f n ). 

3 Examples and numerical results 

In Table 1 we report one example of a second order asymptotic preserving scheme [5]- The scheme is 
also positivity preserving for A < 1 and asymptotically accurate. We also report in Table 2 a third order 
scheme globally stiffly accurate [3] . 

The schemes can be schematically summarized using a double Butcher tableau of the type [llj 



c 


A 




~T 




W 



c 


A 




T 




W 



Note that although the schemes use several implicit evaluations they are still optimal in terms of number 
of evaluation of the collision operator. This, in fact, is characterized only by the number of explicit 
function evaluations. We used the notation name(fc, <te, ctj) where k is the order and o~e, o~i characterize 
the number of evaluations of the explicit and implicit schemes respectively. 

The numerical test is an homogeneous relaxation problem in the two dimensional velocity space. The 
molecules are Maxwellian and a fast spectral method [TU] is used to compute the collision operator with 
N v = 64 grid points in each velocity direction and a grid [— v max , u ma x] 2 with w max = 37r. In this case, 
the exact solution is given by 



f(v,t) 



1 

2irS 2 a 2 



25-1 



l-S v 2 
2S a 5 



exp 



2S<j< 



S(t) = 1 



exp (~cr 2 t/8) 



where we took u = 1. The figure [T] shows the error of the schemes for different choices of the time step 
At (stability condition for the explicit Euler scheme is At = 1). We can clearly observe the expected 
accuracy of the schemes even for large time steps. 
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Table 1: Tableau of the second order IMEX-BE(2,2,4) asymptotic and positivity preserving IMEX scheme. 
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Table 2: Tableau of the third order IMEX-BE(3,5,5) globally stiffly accurate IMEX scheme. 
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